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Static QCD potential at three— loop order 
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We compute the purely gluonic contribution to the static QCD potential at three-loop order. 
This completes the computation of the static potential at this order. 
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For more than 30 years, the static QCD potential 
VQCD(f) has been studied extensively for the purpose 
of elucidating the nature of the interaction between 
heavy quark and antiquark. Generally, Vqcd(?') at 
short-distances can be computed accurately by pertur- 
bative QCD. On the other hand, the potential shape at 
long-distances should be determined by non-perturbative 
methods, such as lattice simulations or phenomenologi- 
cal potential-model analyses or computations based on 
string- inspired models. 

Computations of Vqcd('') in perturbative QCD has a 
long history. At tree- level, Vqcd [r) is merely a Coulomb 
potential, —Cpcts/''' {Cf — 4/3 is a color factor), aris- 
ing from one-gluon-exchange diagram. The 1-loop cor- 
rections (with massless and/or massive internal quarks) 
were computed in P, @|- The 2-loop correction (with 
massless internal quarks) was computed in The 2- 
loop correction due to massive internal quarks was com- 
puted in [1] (partly corrected in ^]).[28)] The logarith- 
mic correction at 3-loop originating from the ultrasoft 
scale was first pointed out in [l| and computed in [7|. 
Renormalization-group (RG) improvement of Vqcd {f) at 
next-to-next-to-leading logarithmic order was performed 
in Q. A logarithmic contribution at 0{a\) was com- 
puted in 0. The contributions of the massless quark 
loops to the 3-loop correction were computed in • The 
only remaining correction at 3-loop order is the purely 
gluonic contribution, which we compute in this paper. 

For a long time, the perturbative QCD predictions of 
VQCDlf) were not successful in the intermediate distance 
region, relevant to the bottomonium and charmonium 
states. In fact, the perturbative series turned out to be 
poorly convergent at r ^ 0.1 fm; uncertainty of the se- 
ries is so large that one could hardly obtain meaningful 
prediction in this distance region. Even if one tries to 
improve the perturbation series by certain resummation 
prescriptions (such as RG improvement), scheme depen- 
dence of the results turns out to be very large; hence, one 
can neither obtain accurate prediction of the potential in 
this region. It was later pointed out that the large uncer- 
tainty of the perturbative prediction can be understood 
as caused by the 0(Aqcd) infrared (IR) renormalon con- 
tained in Vqcd(?') [111] . 

The situation has changed dramatically since the dis- 



covery of the cancellation of ©(Aqcd) renormalons in the 
total energy of a static quark-antiquark pair Etot{r) = 
Vqcd {r) + 2mpoie [13] ■ Convergence of the perturbative 
series for Etotir) improved drastically and much more 
accurate perturbative predictions for the potential shape 
became available. It was understood that a large un- 
certainty originating from the O(Aqcd) renormalon in 
VQCD(f) can be absorbed into twice of the quark pole 
mass 2r7ipoio. Once this is achieved, perturbative uncer- 
tainty of Etot[r) is estimated to be much smaller. 

Then it was readily recognized that perturbative con- 
vergence of Vqcd(?') can be improved by adding a 
(r-independent) constant at each order of the pertur- 
bative expansion, since the ©(Aqcd) renormalon is 
r-independent. The conventional prescription to fix 
VqcdI?") — > at r — > cx) is not optimal as the conver- 
gence of perturbative series is worse at larger r; rather 
it is better to fix Vqcd(?') at some small distance. As it 
turned out, VqcdI?") becomes steeper at r ^ 0.1 fm as 
the order of the expansion is raised, hence convergence 
of perturbative series becomes worse if we fix Vqcd {t) at 
r — > oo. This feature, that the perturbative potential be- 
comes steeper than the Coulomb potential as r increases, 
is understood, within perturbative QCD, as an effect of 
the running of the strong coupling constant [T3|. In fact, 
several studies have shown that perturbative predictions 
for Vqcd {f) agree well with phonomenological potentials 
and lattice calculations of VqcdI?") in the intermediate 
distance region i, [H, [11 [H . 

The improvement of the situation opened up vast ap- 
plications of the QCD potential in heavy quarkonium 
physics [l^ . For instance, higher-order computations 
of VqcdC?") play crucial roles in precise determinations 
of mc,TOb,mt from the masses of charmonium, bot- 
tomonium and (would-be) toponium states. The 3-loop 
correction to Vqcd(?') is one of the missing parts in 
these computations and also in recent efforts to com- 
plete next-to-next-to-next-to-leading order corrections to 
heavy quark production near threshold at e^e~ colliders 
[TtI . Another application is a precise determination of 
as^ from comparison of the perturbative prediction and 
lattice computations of Vqcd [t) [ll] . 

The static QCD potential is defined from an expecta- 
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tion value of the Wilson loop as 

1 {OlTrPe^sfcd^'^^f' |0} 
VQCD[r) = - lim — In- 



47r 



(0 I Tr 1 1 ) 



(1) 



where q = \q\; C is a rectangular loop of spatial extent 
r and time extent T. The second equality defines the 
F-scheme coupling contant, ay(g), in momentum space. 
We employ dimensional regularization with one temporal 
dimension and d — D — 1^3 — 2e spatial dimensions. 
A prefactor is included such that av{q) is defined to be 
dimensionless; je = 0.5772... denotes the Euler constant. 

In perturbative QCD, ay(q) is calculable in series ex- 
pansion of the strong coupling constant. We denote the 
perturbative evaluation of as 



aViq) = as in) ^ 

71=0 



with 



£ = \og{^i/q). 



(2) 



(3) 



Here, asifJ.) denotes the strong coupling constant renor- 
malized at the renormalization scale /x, defined in the 
modified minimal subtraction (MS) scheme; de- 
notes an ri-th-degree polynomial of £. The RG equation 
of as(/i) is given by 



-as 
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n+l 



(4) 



where /3„ represents the (n-l-l)-loop coefficient of the beta 
function. [231 For n < 2, the only part of the polynomial 
Pn {£) that is not determined by the RG equation is a„ = 
P„(0). For n > 3, includes IR divergences in terms 

of poles of £ and assoicated logarithms, whose coefficients 
are not determined by /?„. At 3-loop order, we have 



Psie) = as + (6a2/3o + 4ai/3i + 2ao/?2)^ 

+ (12ai/3o' + 10ao/3o/3i)^' + SaoPo^i^ (5) 
1 
e 



(6) 



Cf = {N^ - l)/{2Nc) and Ca = denote the eigenval- 
ues of the quadratic Casimir operators for the fundamen- 
tal and adjoint representations, respectively, of the color 
SU{Nc) gauge group; TV^ = 3 in QCD. 

The IR divergence is an artifact of the strict pertur- 
bative expansion of Vqcd(?') in as; beyond naive per- 
turbation theory, this IR divergence is absent and reg- 
ularized by the energy difference between color-singlet 
and octet intermediate states. The difference between 
Vqcd(?') and its perturbative expansion [Vqcd('')]pt can 



be treated systematically within the effective field theory 
"potential non-relativistic QCD" [l^. ([Vqcd(?')]pt is 
obtained from ^qcd if we replace ay {q) in eq. H]) by 
ay^{q).) This difference 



Vqcd (f 



us 



PT 



(7) 



is given by contributions of ultra-soft (US) degrees of 
freedom. [3q| In the region r <C Aq^j-,, the leading-order 
contribution to [Vqcd(''')]us in double expansion in as 
and log(a5) is readily obtained from the result of as 



Vqcd(?') 



CfCa^s 



US,LO 
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2 log (Ca as 



+ 81og(^r) 



2 + 67B 



(8) 



Upon Fourier transform, 1/e and log/i terms of eqs. ^ 
and ([8]) cancel each other. In general, at r < 
^QCD' '^^'^ may perform operator-product-expansion 
of [Vqcd('')]us as multi-pole expansion in r. In this 
case, non-perturbative contributions to [Vqcd]us are 
parametrized in the form of non-local gluon condensates. 

We may classify ^3 in ^ according to the powers of 
the number of flavors ni of the internal quarks: 



as 



,3 =(3) 



= (0) 



= (0) 



(9) 



The purpose of this paper is to compute . 

Let us describe our calculational procedure. At tree- 
level and at 1-loop order, computation of [Vqcd('")]pt 
is more or less trivial. The 2-loop correction to 
[^qcd(?')]pt is expressed in terms of 5 master integrals, 
all of which are expressed in terms of F function and ra- 
tional functions of e Hence, we may easily obtain 
expansion coefficients in e necessary for the 3-loop com- 
putation. 

We first generate 3-loop Feynman diagrams for the 
scattering of stati c q uark and antiquark using GRACE 
[131 and QGRAF There are about 20,000 diagrams; 
we confirmed that the diagrams generated by the two 
programs coincide. Next step is to eliminate iterations 
of the lower-order potential at the diagram level, which 
includes appropriate rearrangements of color factors asso- 
ciated with diagrams; we use the general algorithm devel- 
oped in [231] . This procedure eliminates diagrams which 
contain pinch singularities. Subsequently the color factor 
for each diagram is simplified using the program color 
provided in 24]. 

Our computation is carried out in Feynman gauge. 
The loop integrals are classified according to different 
numerators and denominators. At an early stage of the 
computation, we identify those integrals which are triv- 
ially zero in dimensional regularization and eliminate 
them. To reduce the labor of computation, we collect 
integrands with a common denominator and cancel nu- 
merators against denominators as much as possible, by 
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appropriately expressing numerators in combinations of 
factors in the denominator. After these processes, we 
were able to express the 3-loop correction to [Vqcd('^)]pt 
in terms of about 1700 integrals. 

Following the standard procedure of contemporary 
loop computations, these integrals are expressed in terms 
of a small set of integrals (master integrals) through 
the reduction procedure using integration-by-parts (IBP) 
identities [1^. To carry out the reduction efficiently, we 
use the Laporta algorithm [2^. In addition to known 
techniques, we implement some improvement to this re- 
duction algorithm. For instance, we temporarily assign 
a numerical value to D and reduce integrals to sim- 
pler ones using IBP identities. Reduction process com- 
pletes swiftly since manipulation of numerics is consider- 
ably faster than symbolic manipulation involving ratio- 
nal functions of D. We retrace the reduction process and 
identify a minimal set of necessary IBP identities for this 
reduction. Then we reprocess the reduction (without as- 
signing a numerical value to D) using the minimal set 
of identities, after rearranging the order of these iden- 
tities optimally. In the end, the 3-loop correction to 
[^qcd('')]pt is expressed in terms of 40 master integrals. 
All the processes are automatized and the integrals are 
reduced one after another. The reduction processes re- 
quired roughly 3 weeks' CPU time of a comtemporary 
desktop computer with 5 GB memory. 

Out of 40 master integrals, 17 integrals can be ex- 
pressed in terms of F function and rational functions 
of D. The rest of the master integrals are expanded in 
Laurant series in e and their expansion coefficients are 
evaluated analytically if possible and numerically other- 
wise. (For some expansion coefficients, analytical val- 
ues are available in the literature.) Numerical evalua- 
tion of the expansion coefficients are carried out in two 
ways: (a) to evaluate Feynman parameter integrals us- 
ing sector decomposition, and (b) to evaluate integrals in 
Mellin-Barnes representation. Typical relative accuracy 
in numerical evaluation of the expansion coefficients is of 
order 10~^. Details of our computation will be described 
elsewhere. 

Our final result reads 



= (0) 



(502.22(12)) Ci + (- 



136.8(14)) 



^abcd^abcd 

Na 



(10) 



with the color factor d^^'^ 



^abcd , 



/Na = NciN^+6)/48 M- 
For completeness, we combine our result with that of [1(| 
and list the numerical values of 03, defined in eq. (O, for 
TVc = 3 and m = 3, 4, 5 in Tab. D 

At every stage of the computation we performed nu- 
merous cross checks. At every step we have written (at 



ni 


3 


4 


5 


0.3 


5199(3) 


3160(3) 


1460(3) 



TABLE I: Numerical values of as, defined in eq. (|6}, for dif- 
ferent values of nj. 
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FIG. 1: VQCD(r-) = [VQCD(r)]pT + [VQci>{r)]vs up to 0{at) 
and ©(Oglogas): (a) In the toponium region 3 lines, corre- 
sponding to /i = 25, 50 and ICQ GeV, are plotted with ni — 5; 
(b) Comparison with the lattice computations in the quenched 
approximation p3. [27I]. We set n; — 0. The distance region 
corresponds roughly to the size of T(1S) state. 



least) two independent programs and checked that re- 
sults mutually agree. We derived many relations among 
different types of integrals and checked that, when the in- 
tegrals are expressed by master integrals, these relations 
are satisfied. Renormalizability of the QCD potential 
with the known renormalization constant of the strong 
coupling constant, as well as reproduction of known IR 
divergence, provide non-trivial cross checks. We have re- 
produced d'^\ aip , and the coefficient of C|. in a"^^ [l3| 
analytically. We also computed the coefficients of C^, 



CaCf and df'"^df'"^/NA in a^^^ numerically and con- 
firmed that our results agree with those of [loj within the 
estimated errors. The last comparison provides a strong 
cross check on correctness and accuracy of our result (|10p , 
since the expansion coefficients necessary to compute the 
result pO|) are common to the ones necessary to compute 



-(3) -(2 



3 , U.3 , u-3 except for a single coefficient. 
Now we compute Vqcd(?'), as given by the sum of 
[Vqcd(7')]pt and [V'QCD(r)]us, including all the correc- 
tions up to 0{a^) and ©(a^ log 05). Namely we use the 
series expansion ^ wp to n — 3 for the former and eq. ([5]) 
for the latter. 

In Fig.[lja) we plot our prediction for Vqcd(?') in the 
distance region corresponding to (would-be) toponium 
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states. 3 lines are plotted, corresponding to ^ = 25, 50 
and 100 GeV, with = 5 and as{Mz) = 0.1176. We 
added a constant to each prediction such that it takes 
a common value at r = 0.01 GeV~^. The difference of 
the 3 lines are hardly visible, showing stability of the 
prediction. 

In Fig. [IJb) we compare our prediction with the lat- 
tice data in the quenched approximation [iJi [13|- 
cordingly we set m = 0. We used the central value of 
roA|j^°P = 0.574 ± 0.042 [3 to fix the relation between 

the lattice scale and A|j|!°p, where tq denotes the Som- 
mer scale. Hence, the only adjustable parameters in our 
comparison are r-independent constants to be added to 
the potentials, whose values are chosen such that all the 
potentials coincide at rA|j^°^ = 0.1. It is customary 
to interpret ro = 0.5 fm when comparing this scale to 
one of the real world. Roughly the potential shape in 
the displayed range r < ro/2 accounts for formation of 



the T{1S) state. We plot 3 lines with the scale choices 
A|_|°p//^ = 0.14, 0.07 and 0.035. (The corresponding 
values of as(/i) are 0.216, 0.165 and 0.135, respectively.) 
There is a small but visible dependence on the scale. The 
level of agreement with the lattice data shows that our 
prediction of the potential at this order is good enough 
to warrant quantitative description of the nature of the 
T(IS') state. We confirm the observation that, as we in- 
clude higher-order corrections, agreement of the pertur- 
bative prediction and lattice computations improve s up 
to larger distances. According to the analyses in (T3l.ll8j|. 
we anticipate that the agreement would get even better 
if we resum logarithms via RG, or, appropriately choose 
the scale /i as a function of r, provided that the IR renor- 
malon is subtracted. 
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